Schooling behavior driven complexities in a fear-induced prey–predator system with harvesting under deterministic and stochastic environments

The well-being of humans is closely linked to the well-being of species in any ecosystem, but the relationship between humans and nature has changed over time as societies have become more industrialized. In order to ensure the future of our ecosystems, we need to protect our planet’s biodiversity. In this work, a prey–predator model with fear dropping prey’s birth as well as death rates and nonlinear harvesting, is investigated. In addition, we consider that the consumption rate of predators, i.e., the functional response, is dependent on schooling behavior of both species. We have investigated the local stability of the equilibrium points and different types of bifurcations, such as transcritical, saddle-node, Hopf and Bogdanov–Takens (BT). We find that consumption rate of predator, fear and harvesting effort give complex dynamics in the neighbourhood of BT-points. Harvesting effort has both stabilizing and destabilizing effects. There is bistability between coexistence and predator-free equilibrium points in the system. Further, we have studied the deterministic model in fluctuating environment. Simulation results of stochastic system includes time series solutions of one simulation run and corresponding phase portraits. Notably, several simulation runs are conducted to obtain time series solutions, histograms, and stationary distributions. Our findings exhibit that during stochastic processes, model species fluctuate around some average values of the deterministic steady-state for lower environmental disturbances. However, higher values of environmental disturbances lead the species to extinction.

the main consideration in prey-predator models 4,18 . However, in reality, fear of predators can affect more than reproduction; it can also affect the death rates of prey population 5,19,20 . There is no doubt that harvesting is an influential issue from both an ecological and economic standpoint as well as from a social perspective. Moreover, there is a possibility that predators and/or prey species can exploit the ecosystem when there is an abundance of food. As a result, it is necessary to harvest species. There have been a number of researchers who have explored this area extensively in recent years [21][22][23] .
Stochasticity (i.e., variation in the model parameters due to random effects) in the population models improves the accuracy, realism, and utility of the models. Generally, population dynamics are affected by two types of stochasticity 24,25 . They are (i) Demographic stochasticity: random variation in the number of births and deaths in a population caused by the discrete nature of individuals. (ii) Environmental stochasticity: variability on the environmental conditions such as temperature, humidity, pH, rainfall, etc. It is obvious that survival and reproduction tend to be affected by environmental conditions. Environmental stochasticity applies to both small and large populations, whereas Demographic stochasticity negligible in case of large populations. Additionally, it is desirable to measure the variability of outcomes within a conservation or restoration framework 26 . To expose the actual dynamics of a population model within open environment, environmental stochasticity should always be considered as it is impossible to keep environmental conditions constant over time. Incorporation of environmental fluctuations or demographic stochasticity into the modeling approach are important components. In several existing literature 27 , author has shown that continuous fluctuation of environmental conditions can lead to random fluctuation in the important model parameters to a greater or lesser extent. They are mainly birth rates, death rates, carrying capacity, competition coefficients and all other parameters involved in the dynamical system. It has been established that environmental noise has a significant effect on deterministic systems when random disturbances are introduced 22,28-31 . Noise from the environment adversely affects almost all ecosystems. Thus, prey-predator models cannot ignore the shifting environmental effects, whereas stochastic models can accurately predict dynamics when the environment changes 22 . Belabbas et al. 32 investigated a new approach of a stochastic prey-predator model with protection zone for the prey and found rich dynamics of the system.
Following the above discussions, we were motivated to visit the state as discussed here. There are hardly any studies considering predator-dependent functional response describing both predatory and prey schooling behaviors 10,33 . Additionally, a limited number of literatures address the issue of fear effect affecting the death of prey population 5,34,35 . Here, we intend to explore the deep insights of a prey-predator model with Cosner-type functional response 10 , fear that affects the growth and death of prey population, and nonlinear harvesting of the predator population. To the best of our knowledge, none yet studied the combined effects of double fear, nonlinear harvesting on Cosner-type functional response to fill the gap in extant research. Additionally, we incorporate white noise 36,37 due to the perturbation of environmental conditions. Our second objective is to determine how environmental noise affects the dynamics of the system. Moreover, a numerical comparison between deterministic and stochastic models is made.

Deterministic model
In a region under consideration, let at any instant t > 0 , x and y represent the prey and predator population densities, respectively. The rate of change of each model species density at time t is made on the following assumptions: 1. Prey population grow logistically in the absence of predator with birth rate r, which is affected by the fear ( f 1 ) of predator (when predators are around). 2. There is a reduction in the rate of prey density change due to three types of death, namely, natural death with the rate d 1 , fear related death 5 with the level of fear f 2 and over crowding death with the rate d 2 . 3. Also, the rate of change of prey density decreases due to predation of predator population following a predator-dependent functional response describing both predatory and prey schooling behaviors 10 . Response function is expressed in functional form describing as ζ(x, y) = cxy 1+chxy , where c denotes the rate of consumption and h represents handling time of predator for one prey. 4. Predator population survive in the system by consuming prey population only. They grow with conversion efficiency c 1 of prey biomass into predator biomass. 5. Predator population harvested from the system which reduces its rate of density. We consider a nonlinear harvesting term (Michaelis-Menten type) given by, H(y) = qEy p 1 E + p 2 y . Here, parameters q and E, respectively, represent the catchability rate and harvesting effort. It is easy to observe that H → q p 1 y as E → ∞ for a fixed value of y. Also, H → q p 2 E as y → ∞ for a fixed value of E. Therefore, at higher effort levels, p 1 is proportional to the stock level-catch rate ratio and at higher levels of stock, p 2 is proportional to the effort level-catch rate ratio. 6. Lastly, we assume that the predator population experience natural as well as over crowding related death with the rates d 3 and d 4 , respectively.
Keeping all these above assumptions in mind, we formulate the following prey-predator model: (1) www.nature.com/scientificreports/ System (1) is to be analyzed with the initial conditions x(0), y(0) > 0 . All the model parameters are assumed to be positive constants and their hypothetical values that we used for numerical calculations are as follows: In Table 1, we have provided system's equilibria, sufficient conditions of their existence and stability. Mathematically, it is difficult to determine the existence of coexistence (interior) equilibrium point(s) by the given nullclines. So, we visualize it numerically (see Fig. 1). It is apparent from the figure that on increasing the value of E, number of coexistence equilibrium points reduces and after a certain range there is no coexistence equilibrium point. Table 1, it is clear that the equilibrium E 0 is stable if r < d 1 , which is opposite to the existence condition of E 1 . That is, equilibrium E 0 is stable whenever the equilibrium E 1 does not exist, and hence these two equilibria are related via transcritical bifurcation. Using Sotomayor theorem 38 , we can easily prove that model system (1) experiences transcritical bifurcation at the trivial equilibrium point E 0 as the growth rate of prey crosses a critical value r [TB] = d 1 .

Transcritical bifurcation. From
We visualize the transcritical bifurcation graphically in Fig. 2. It is clear from the figure that when the value of r is less than d 1 = 0.1 , the equilibrium point E 0 only exists and it is stable. If we increase the value of r ( r > r TB = d 1 = 0.1 ) then E 0 becomes unstable and the equilibrium point E 1 exists and becomes stable. Table 1. Sufficient conditions for the existence and stability of different equilibrium points of system (1).

Equilibria
Condition of existence Stability condition x * and y * are the positive solution(s) of the nullclines   www.nature.com/scientificreports/ Hopf bifurcation. One of the most common dynamics in interacting population dynamics is oscillating behavior, which implies that there is a Hopf bifurcation. By local changes in equilibrium properties, Hopf bifurcation describes when a periodic solution appears or disappears. In this section, we study the Hopf bifurcation through the coexistence equilibrium E * with respect to the model parameter E. Discussion for the existence of Hopf bifurcation is as follows: As it is easy to follow, we verify Hopf bifurcation numerically. We have considered the parameters value same as in (2) Therefore, the transversality condition for Hopf bifurcation is also satisfied at E = E [HB] . Thus, these results confirm that the system (1) experiences a Hopf bifurcation 2 around E * (2.618402886, 2.352228027).
Saddle-node bifurcation. Sotomayor's theorem says that saddle-node bifurcation may occur at coincident equilibrium points depending on threshold value of the bifurcation parameters. Since the analytical result is difficult to follow, we verify the existence of saddle-node bifurcation of system (1) numerically for given set of parameters value (2)

Bogdanov-Takens bifurcation. A number of one-dimensional bifurcations such as transcritical, Hopf
and saddle-node have been studied in earlier segments. Each of these bifurcations belongs to one parametric bifurcation. In the present segment, we will discuss two parametric bifurcation, i.e., of codimension two bifurcation. Bogdanov-Takens bifurcation (BT) is a such type of bifurcation. When Hopf and saddle-node bifurcation curves collide, a bifurcation of this type occurs at the vicinity of colliding point. In this case, the critical values of two bifurcation parameters give zero eigenvalues of multiplicity two in the Jacobian matrix of the system.
One can used Kuznetsov's 39 technique to achieve the standard form of BT-bifurcation. Due to model complexity, we did not find any explicit expression of BT-bifurcation. So, we verify it numerically, which is discussed below.
System (1) undergoes BT-bifurcation for the set of parameters given in (2) except c and E. We consider c and E as bifurcation parameters and observe that saddle-node and Hopf bifurcation curves collide with each other at (E, , and the instantaneous equilibrium points are (11.12253608, 1.485272279) and (13.778434, 1.13304881), respectively. Therefore, system (1)  Numerical results of the deterministic system (1). To explore the rich dynamics of the system (1), we perform extensive numerical simulations in this section. For numerical simulations, we choose a set of biologically feasible hypothetical parameter values as given in (2).
A system may exhibit rich dynamical behavior if its dynamics are investigated in different bi-parameter spaces. So, here we plot two bi-parametric bifurcations (Figs. 3a, 5a). In Fig. 3a, we plot saddle-node, Hopf and Bogdanov-Takens bifurcations in E − c plane, which divide the entire region into six different regions ( R 1 − R 6 ) of different dynamical behaviors. Complete phase portraits are drawn for every region in order to understand their dynamics transparently, Fig. 3b-i. It is apparent from Fig. 3a that the saddle-node curve divides the whole region into two parts. Among these two parts one part is region R 1 , which has no interior equilibrium point (see Fig. 3b) and other part is sum of rest regions ( R 2 − R 6 ), that contain two interior equilibrium points (see Fig. 3c-i). Next, for lower and higher values of E, there exist Hopf bifurcation curves, which divide two portions of the region into two subregions that are ( R 2 /R 3 ) and ( R 5 /R 6 ). In one side of the Hopf curve, one of the two interiors is stable spiral and other one is saddle, Fig. 3c,i (correspond to the regions R 2 & R 6 , respectively). On the other hand, as the value of E and c crosses Hopf curve, stable spiral point becomes unstable spiral (stable limit cycle) and saddle one remain same, Fig. 3d,h (correspond to the regions R 3 & R 5 , respectively). Similar as Hopf curve, there exist two Homoclinic curves for lower and higher values of E, which divide a subregion into two another subregions ( R 3 /R 4 ) and ( R 4 /R 5 ). As the values of (E, c) closest to the Homoclinic curve system gives larger amplitude-limit cycle and on the curve it gives maximum amplitude-limit cycle, Fig. 3e,g. Further, if the values of (E, c) crosses the Homoclinic curve, stable limit cycle destroys and the equilibrium point becomes unstable (see Fig. 3f corresponds to the region R 4 ). Also, it is clear from Fig. 3a that there exist two BT-points ( BT 1 = (1.428292550, 0.1054627152) www.nature.com/scientificreports/ and BT 2 = (7.345497850, 0.2137076372) ), where aforementioned three bifurcation curves meet with each other. Therefore, the neighbourhood of BT-points exhibits complex dynamics of the system. It is observed that the system preserves bistability for the regions R 2 , R 3 , R 5 and R 6 . For a fixed set of parameter values, depending on the initial population size the trajectories initiating inside the green regions converge to the stable interior (stable limit cycle) whereas the trajectories initiating in the red regions converge to the predator-free equilibrium point, Fig. 4. In Fig. 5a, we plot another bi-parametric bifurcation in f 1 − q plane, which is divided into four subregions ( R 1 − R 4 ) of different dynamics by different bifurcation curves. Phase portraits of every region are plotted (see Fig. 5b-f). We observe that there is no interior for diagonally higher values of both of the parameters f 1 and q whereas two interiors exist when the values of f 1 and/or q are lower, Fig. 5a. Detailed discussions of the figure are same as previous one (Fig. 3).
Next, to observe the effect of important model parameters explicitly, we plot one parameter bifurcation diagrams while other parameters are fixed, Fig. 6. It is clear form Fig. 6a,d,e that for lower values of the parameters r, c and c 1 , there is no interior. As the value of these parameters surpasses saddle-node bifurcation point there exist two interior one of which always remain saddle and another one switches its stability (stable spiral ⇒ unstable spiral) through Hopf bifurcation point. Therefore, up to a certain range of these parameters have destabilizing effect. On the other hand, for lower values of the parameters f 1 , f 2 and E there exist two interior one of them is always saddle and another one switches its stability (unstable spiral ⇒ stable spiral) through Hopf bifurcation  www.nature.com/scientificreports/ point (see Fig. 6b,c,f). But, when the value of these parameters crosses saddle-node bifurcation point, no interior equilibrium point exists. Note that up to a certain range of these parameters, f 1 and f 2 have stabilizing effect whereas E has both destabilizing as well as stabilizing effects.
Interpretation of the deterministic results in the context of biology. For the deterministic system, we mainly investigate different types of bifurcation results. They are Transcritical, Hopf, Saddle-node and Bogdanov-Takens. In the context of a biological system, a bifurcation occurs when a small smooth change made to the parameter values (the bifurcation parameters) of a system causes a sudden ' qualitative' or topological change in the behavior of the system. Bifurcations mainly describe changes in the stability and/or existence of fixed points (equilibrium points) and the bifurcation parameters act as a control parameter. Changes in the control parameter eventually changed the qualitative behavior of the system. Above mentioned bifurcations of the considered system (1) were found for the model parameters r, f 1 , f 2 , c, c 1 , q and E. Therefore, all these ecological parameters (i.e., growth rate of prey, fear, capture rate, harvesting effort) act as control parameters of the proposed system. In Transcritical bifurcation, there is a critical value of the bifurcation parameter from which two equilibrium points exchange their stability. Moreover, on the one side of the critical value one of the two equilibrium points  www.nature.com/scientificreports/ exist while on the other side both equilibrium points exist. Therefore, we have a critical value of an ecological parameter from which species persistency and stability can be described. In case of Hopf bifurcation, there is a critical value of a parameter from which densities of model species either fluctuate from fixed stable densities or vice versa. In the case of Saddle-node bifurcation, coexistence equilibrium points exist or diminish in pair by varying bifurcation parameter i.e., from the critical value of the bifurcation parameter coexistence equilibrium points exists or diminish in pair. Lastly, BT-bifurcation, the point where saddle-node curve and Hopf curve meet is known as BT-point. Thus, all these phenomena occur around the BT-point, i.e., near the BT-point complex dynamics (like species persistency, stability, extinction) occurs. Thus, studying the bifurcation in prey-predator systems has great importance from environmental perspective. In order to maintain an ecological balance between prey and predator populations, identification of bifurcation parameters plays a crucial role in determining an effective control strategy.

Stochastic model
There is no doubt that the system (1) is derived based on the assumption that all the input variables follow deterministic laws and are deterministic functions of time. However, in mathematical modeling of ecosystems, the deterministic system has its limitations since it cannot capture the influence of random environmental fluctuations in its parameters 31 . In order to study the effects of environmental fluctuations on the entire ecosystem, it is reasonable to introduce the noise term into the deterministic model. Mao et al. 40 demonstrated that one or more system parameter(s) can be perturbed stochastically with white noise term to derive stochastic system. Note that the approach to formulate the stochastic model based upon existing deterministic model is not unique 41,42 .
Introducing multiplicative noise terms into the growth equations of both prey and predator populations, we formulate the following stochastic model.
where σ i (i = 1, 2) represent the intensity of environmental fluctuations and B i (t) are standard Brownian motions. Throughout the analysis, we take (�, F , P) as a complete probability space with a filtration {F t } t∈R + satisfying the conventional condition, namely right continuity and increasing, whereas F 0 consists of all P − void sets 40 . Any solution of system (3) subjected to the positive initial condition is an Itô process 43 . Without any loss of generality we assume σ 1 , σ 2 > 0.
Numerical results of stochastic system (3). In this section, we simulate the stochastic model (3) to explore the dynamics of different noise intensities.
In Fig. 7, we plot deterministic and stochastic time series solutions and their corresponding phase portraits. Random perturbations in the model parameters (i.e., environmental noise) destroy the stability of deterministic equilibrium and lead to weak stochastic stability known as stationary distributions. According to Fig. 7a, in the absence of environmental noise, model variables approach to their equilibrium values, and in the presence of environmental noise, stationary distributions are observed. Our results indicate that when noise intensity is relatively low, the stochastic system still maintains some stability. Now, if we increase the value of σ 1 to 0.6 from 0.01, then after some initial transients predator population extinct. However, it is important to note that in this case the prey species is not fluctuating around the deterministic steady-state rather it fluctuate above its equilibrium value (obviously, it happen due to predator extinction) (see Fig. 7b). Next, we choose σ 1 = 2.5 and σ 2 = .01 , in this case both the model species go to extinction, Fig. 7c. Moreover, it is clear from the figure that prey species goes to extinction first, before predators. However, the model species approach to their equilibrium values in the absence of environmental noise. It is important to note that in all of these instances of species extinction, extinction time may vary from simulation to simulation, but extinction is confirmed in every simulation. Lastly, we fixed our deterministic system in oscillatory state and see the effect of environmental disturbances. We observe that lower strength of environmental noise propel the species to fluctuate around the limit cycle (see Fig. 7d). However, higher strength of disturbances (i.e., environmental noise) ultimately lead predator population to extinction after some initial transient dynamics (see Fig. 7e). In this case, fluctuations of prey density occur above the deterministic density.
To be more transparent of the effect of stochasticity, 200 simulations are plotted in Fig. 8 as time series solutions. It can be easily observed that all the solution trajectories fluctuate around the deterministic steady-state due to environmental disturbances, Fig. 8b. Histogram plots of 1000 simulation runs also show that these fluctuations are observed at the stationary distributions, where prey population is distributed within the range (9, 11) and predator population within the range (1.75, 2) (see Fig. 8c).
Population fluctuation due to environmental noise is also reflected in the stationary distributions. Therefore, for better presentation and understanding the stochastic effect on population dynamics, stationary distributions of different noise intensities are attained. The result of stationary distributions obtained from 500 simulations at t = 100 , which is presented in Fig. 9. We observe that the population distributions occur in wider range as the noise intensity increases. Whenever the strength of environmental disturbances is low, prey population are distributed within (2.6, 3.1) and it is estimated that the predator population are within the range of (2.35, 2.55). In the next case they are within (0, 8) and (1.5, 4), respectively for moderate strength of noise. Lastly, for higher strength they are within (0, 12) and (1,5), respectively. In these cases, there are no extinction scenarios associated with these parameters, since they do not satisfy the conditions of extinction. Thus, we conclude from these   www.nature.com/scientificreports/ in most countries is always varying. . . . " Therefore, stochastic models are preferred over deterministic ones. From numerical results we observe that if the intensity of environmental noise is small, prey and predator populations will survive for a long time. Therefore, intensities of environmental noise are conducive to the survival of the population. Also, it is evident that there is an important role of environmental noises for the persistence of prey and predator populations weekly as well as strongly. We notice that the prey population can persist for low intensity of noise, higher growth rate of prey, lower level of both fears and consumption rate. The predator population can persist for higher values of conversion efficiency of predator and lower values of catching capability and low intensity of noise on predator population. Population extinction is a serious issue from an ecological perspective. It is observed that high intensity of noise plays an crucial role in extinction of the species. The extinction of prey species drive predator species towards extinction as they depend on prey species only. Further, under a circumstances when the environmental disturbance for prey is low while for predator is high, only predator population extinct. Furthermore, at the simulation time, we notice that increasing intensity of noise decreases average extinction time. More importantly, chaotic variation in the population size due to environmental fluctuations have serious implication, it reduces species extinction 45 . These findings are consistent with the reality which imply that in any natural ecosystem, environmental fluctuations have significant influences on the persistence and extinction of interacting species. As a result, a conservationist may be able to take a variety of measures to prevent species extinctions by identifying different routes to extinction.

Discussion and conclusion
We propose a prey-predator model where schooling behaviour of both prey and predator are considered by a predator-dependent functional response. Predator-induced fear is assumed to affect the birth as well as death rates of the prey population. In order to conserve resources and manage the social environment, predator population are harvested. The complete investigation or discussion of the system is mainly devoted to the important ecological factors, namely, the fear factors, consumption rate of predator (which depend on schooling behavior of the species) and harvesting. Firstly, we study the dynamics of the deterministic system and then stochastic system. Both these systems were studied in detail. For deterministic system, we have established the existence and stability conditions of the equilibrium points. Different types of bifurcations are also numerically investigated, including transcritical, saddle-node, Hopf, and Bogdanov-Takens bifurcations. We have plotted two parameter bifurcation diagrams in which the aforementioned bifurcation curves are plotted which divide the whole region into subregions of different dynamics. Complex dynamics are observed around the BT-points for harvesting effort, fear, and predator consumption rate. We find that growth rate of prey, consumption rate of predator and conversion efficiency of prey biomass into predator biomass have destabilizing effects. In contrast, fear affecting growth and death of prey has stabilizing effects. The harvesting effort, however, has both stabilizing and destabilizing effects. Further, for lower values of predator consumption and higher values of harvesting effort, there is no interior (coexistence) equilibrium point. A BT-point has been observed around a low harvesting effort and moderate consumption rate, exhibiting complex dynamical behavior around it (BT-point). Therefore, these parameter values compensate with each other to exhibit complex biological dynamics. Moreover, moderate values of both fear affecting prey's growth and catching capability enrich the system with rich biological dynamics around the BT-point. The predator population is at risk of extinction due to higher harvesting effort and lower consumption rates. Higher levels of both fear and catching capability also eliminate predator populations from ecosystems. Both species, however, coexist and switch between stability for lower and higher levels of harvesting effort. Prey and predator populations always live in harmony at lower values of catching capability, regardless of the value of fear dropping prey's birth. As the amount of fear increases, coexisting species change their stability from unstable spiral to stable spiral. Furthermore, we incorporate multiplicative noise terms in the deterministic system to understand the dynamics in the presence of environmental driving forces. We can obtain that there exists a unique positive global solution for the stochastic model and determined the conditions under which a species is likely to persist or fail. We did not include formal mathematical analysis of these results. According to our numerical findings, species survival is closely linked to the intensity of environmental fluctuations. Recently, Rogers et al. 46 showed that "Chaos is not rare in natural ecosystems". Therefore, chaotic nature in the population densities due to environmental fluctuations is an important result. Biologically, chaotic variation in the population size have serious implication: it has the ability to reduce species extinction 45 . Moreover, we can obtain parametric conditions mathematically under which stationary distributions exist in the stochastic system. According to simulation results, these parametric restrictions will not hold for large-amplitude environmental noise. Consequently, it can destabilize the system, and, in that case, no stationary distribution can be found 32,47 . As a result of the existence of stationary distributions, there can be some degree of stochastic stability. In terms of biology, this indicates that both prey and predator populations coexist on a long-term basis, leading to the conclusion that the system is permanent. Based on numerical results of the stochastic system, we have demonstrated that for lower levels of noise, stationary distributions are attained, while high levels of noise lead to species extinction. It is possible to observe that there are two different scenarios of extinction: the first case is the both populations extinct; second case is the only prey population persist while predator extinct. From an ecological viewpoint, comparing the stochastic results with corresponding deterministic result, we observe two interesting facts. For low intensity of environmental fluctuations both prey and predator populations coexist in the long run, i.e., the system is permanent. Another fact is that none of the important ecological factor (growth rate of prey, fear, harvesting, consumption rate of predator) can avoid the extinction of model species when the nature exhibits large amount of environmental fluctuations. Although, all these important ecological factors have significant impact on species persistence and extinction in constant environment (i.e., deterministic environment). Following articles [48][49][50] , we can also define stochastic Hopf bifurcation 48 www.nature.com/scientificreports/ in the considered system. Study of bifurcations in the stochastic framework is an interesting topic. In the near future, we will discuss this topic in more detail.
Human or animal involvement. The current study did not involve any animal or human experiments.